Nonlinear response of dense colloidal suspensions under oscillatory shear: 
Mode-coupling theory and FT-rheology experiments 
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Using a combination of theory, experiment and simulation we investigate the nonlinear response 
of dense colloidal suspensions to large amplitude oscillatory shear flow. The time-dependent stress 
response is calculated using a recently developed schematic mode-coupling-type theory describing 
colloidal suspensions under externally applied flow. For finite strain amplitudes the theory generates 
a nonlinear response, characterized by significant higher harmonic contributions. An important fea- 
ture of the theory is the prediction of an ideal glass transition at sufficiently strong coupling, which 
is accompanied by the discontinuous appearance of a dynamic yield stress. For the oscillatory shear 
flow under consideration we find that the yield stress plays an important role in determining the 
non linearity of the time-dependent stress response. Our theoretical findings are strongly supported 
by both large amplitude oscillatory (LAOS) experiments (with FT-rheology analysis) on suspen- 
sions of thermosensitive core-shell particles dispersed in water and Brownian dynamics simulations 
performed on a two-dimensional binary hard-disc mixture. In particular, theory predicts nontrivial 
values of the exponents governing the final decay of the storage and loss moduli as a function of 
strain amplitude which are in excellent agreement with both simulation and experiment. A consis- 
tent set of parameters in the presented schematic model achieves to jointly describe linear moduli, 
nonlinear flow curves and large amplitude oscillatory spectroscopy. 
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I. INTRODUCTION 

A standard method to probe the viscoelastic charac- 
ter of a material is to measure the time dependent stress 
response to an externally applied oscillatory shear field 
PQ. The simplicity of oscillatory shearing experiments 
presents distinct practical advantages when compared to 
other flow protocols and thus makes desirable a system- 
atic method for the rheological characterization of a ma- 
terial on the basis of the periodic stress response alone. 
For small strain amplitudes the shear stress is a simple 
harmonic function, oscillating with the fundamental fre- 
quency dictated by the applied strain field. The details of 
the microscopic interactions underlying the macroscopic 
stress response are encoded in the familiar storage (C) 
and loss (G") moduli of linear response. General aspects 
of the viscoelastic character of the material can thus be 
inferred from the magnitudes of the moduli as a function 
of frequency. 

While, for many systems of interest, the linear response 
regime is well understood, for practical applications, such 
as the production and processing of materials in industry 
[2], it is necessary to consider deformations of finite, of- 
ten large, amplitude. In the nonlinear regime, the stress 
response to a sinusoidal excitation contains higher har- 
monic contributions, which arise from the nonlinearity of 
the underlying constitutive relation expressing the stress 
as a functional of the strain [3H5] . For many complex ma- 
terials, consideration of the fundamental frequency alone 



proves insufficient for describing the physical mechanisms 
at work for finite strain amplitude. Analysis based purely 
on the linear complex modulus as a function of frequency 
can thus be expected to give only a partial mechani- 
cal characterization of the system under study (see e.g. 
[6l[8]). This failing is found to be particularly pronounced 
for yield stress materials such as aqueous foams [9] and, 
as we will argue in the present work, colloidal suspen- 
sions close to, or beyond, the point of dynamical arrest. 
Although such systems are predominately elastic in char- 
acter, they exhibit a complex transient response to oscil- 
latory shear in which the viscous dissipation mechanism 
present at small strain amplitudes crosses over to a plas- 
tic flow as the amplitude is increased. The nonlinear 
stress response reflecting the onset of plastic flow gives 
rise to a strong increase in the amplitudes of the higher 
harmonics. 

The emerging discipline of Fourier transform rheol- 
ogy (FT-rheology), originating in the work of Wilhelm 
and co-workers (see, e.g. [4j [101 Q2] ) , aims to quan- 
tify the nonlinear response of complex fluids by analyz- 
ing the harmonic structure of the stress signal measured 
in large amplitude oscillatory shear (LAOS) experiments 
(for recent developments see [13). Despite considerable 
progress on the experimental side, the theoretical descrip- 
tion of the nonlinear regime remains unsatisfactory. The- 
oretical treatments capable of capturing higher harmonic 
contributions have been largely restricted to phenomeno- 
logical models based on the ideas of continuum rheology 
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[3 H HOI HD [H [15]. A more refined description of the 
nonlinear response is provided by mesoscopic models in 
which the time evolution of explicit coarse-grained de- 
grees of freedom is governed by specified dynamical rules 
[ToT - fT8] . While such approaches are capable of capturing 
generic features of the response, they are not material 
specific and make no explicit reference to the underlying 
particle interactions. 

Recently, progress in making the connection between 
microscopic and macroscopic levels of description has 
been made for the case of dense colloidal suspensions 
subject to time-dependent flow [T9l 120] . The develop- 
ments in classical nonequilibrium statistical mechanics 
presented in [115] and [20] extend earlier work focused on 
the simpler, but fundamental, case of steady shear flow 
[23 [23]. The mode-coupling-type approximations em- 
ployed in [T3 [23 [23 123] capture the slow structural re- 
laxation leading to dynamical arrest in strongly coupled 
systems (i.e. dispersions at high volume fraction or with 
a strongly attractive potential interaction), with the con- 
sequence that the macroscopic flow curves 17(7) attain a 
finite value in the limit of vanishing rate, for states which 
would be glasses or gels in the absence of flow. The fi- 
nite value of the stress in the slow flow limit identifies 
the dynamic yield stress. The relationship between the 
dynamic yield stress and its more familiar static coun- 
terpart is analogous to that between stick and slip fric- 
tion in engineering applications. A prediction of particu- 
lar importance made by the the MCT-based approaches 
of [13 [23 123 123] is that the dynamic yield stress ap- 
pears discontinuously as a function of coupling strength, 
in clear contrast to mesoscopic models [13 HZ] which pre- 
dict a continuous power law dependence. The notion of 
yield stress was considered in a more general and ab- 
stract sense in [24 , in which a dynamic yield stress sur- 
face, describing yielding under more general non-shear 
deformations, was calculated (see also [25]). 

Although the closed, microscopic constitutive equa- 
tion presented in [20] is of considerable generality, the 
combined difficulties of a large time-scale separation be- 
tween microscopic and structural relaxation times, spa- 
tial anisotropy, and lack of time-translational invari- 
ance presented by many problems of interest makes di- 
rect numerical solution of the equations impossible at 
the present time. In order to both facilitate numeri- 
cal calculations and expose more transparently the es- 
sential physics captured by the fully microscopic theory 
of [20 a simplified 'schematic' model has been proposed 
[24j . Schematic models have proved invaluable in the 
analysis and assessment of microscopic mode-coupling 
approaches, both for quiescent systems [2S| and under 
steady shear flow |27j . in each case providing a simpler 
set of equations which aim to retain the essential mathe- 
matical structure of the fully microscopic theory. While 
the schematic model reduction performed in [24] leads 
to loss of the 'first principles' character of the approach, 
the mathematical connections between full and schematic 
theory nevertheless serve to elevate the schematic model 



above purely phenomenological approaches. 

In the present work we will consider application of 
the schematic model derived in 24 to the problem of 
large amplitude oscillatory shear. Although the tensorial 
schematic model of [23 is closely related to the earlier 
F^ 2 model derived in [27], application of the tensorial 
model to a simple shear flow geometry does not exactly 
reproduce the F± 2 model. The study of time-dependent 
flows, not considered in earlier work, revealed that cor- 
rections to the original F^ 2 model were necessary to cap- 
ture correctly the response to rapidly varying flows. The 
modifications thus introduced lead to small differences 
in the steady state rheological predictions. Nevertheless, 
the present schematic models describes the same phe- 
nomenology as the previous model |28) when applied to 
steady shear. 

Comparison of theoretical predictions with experimen- 
tal data for thermosensitive core-shell particles, dispersed 
in water, has been performed using the F^ 2 model [2"Tj . 
These particles have the very convenient feature that the 
volume fraction of the system may be varied continuously 
over a considerable range, simply by tuning the temper- 
ature of the system. Moreover, the finite polydispersity 
in particle size effectively suppresses crystallization, such 
that studies of dense fluid and glassy states are not com- 
plicated by an intervening fluid-crystal transition. In a 
series of works, theory and experiment have been com- 
pared for the flow curves under steady shear [23130] and, 
more recently, for both flow curves and linear response 
moduli [3U [32] ■ A particular strength of the F^ 2 model 
(inherited by the more recent model of [21]) is that both 
flow curves and linear viscoelastic moduli can be simulta- 
neously and accurately fitted over many decades of shear 
rate and frequency, respectively, using a consistent and 
physically meaningful set of fit parameters. In |32j a 
combination of experimental techniques were employed, 
which enabled measurement of the flow curves and linear 
response moduli over eight and nine orders of magnitude 
in shear rate and frequency, respectively |32j . Although 
certain discrepancies between experiment and theory at 
low frequencies remain to be fully understood, the gen- 
eral level of agreement is impressive. Reassuringly for the 
schematic models, the complete microscopic MCT calcu- 
lations possible for the linear response moduli agree with 
the data from the monodisperse samples on the 40% error 
level [3PJ. 

The nonlinear rheology of thermosensitive microgel 
particles (identical to those considered in the present 
work) was addressed in a recent experimental study, fo- 
cused on the stress response to steady and large ampli- 
tude oscillatory flow [31]. In addition to a study of the 
stress overshoot following the onset of shear flow (see 
also [3SJ), both the strain dependence of the storage and 
loss moduli and the higher harmonic contributions were 
analyzed. Despite employing the same thermosensitive 
particles and LAOS flow protocol, the study [31] should 
be regarded as complementary to the present work. In 
[34] volume fractions well above random close packing 
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were investigated (0 > 0.64), suggesting considerable de- 
formation of the particles themselves, whereas we focus 
here on packing fractions around the glass transition. 
Moreover, emphasis in the present work is placed on as- 
sessing the MCT based schematic theory presented in 
[24] for a nontrivial flow history, namely large amplitude 
oscillatory shear, and comparison of the theoretical pre- 
dictions with experiment. This comparison provides the 
first truly time-dependent test of this recently developed 
schematic model beyond the simple case of step strain 
already considered in [2"4"j . 

The paper will be organized as follows: In Section |TT] 
we summarize the microscopic starting points underly- 
ing our theoretical approach, before proceeding to give a 
compact overview of the linear and nonlinear response of 
viscoelastic systems, relevant for the subsequent analysis. 
In Section IIIII we introduce the schematic MCT model 
and discuss its relation to previous work. In Section |IV| 
we discuss the Brownian dynamics simulation algorithm 
used to generate results supplementary to those of theory 
and experiment. Section [V] contains the experimental de- 
tails. In Section [VI] we first present purely theoretical re- 
sults, in order to establish the phenomenology predicted 
by the schematic model. We then consider the results 
of our two-dimensional simulations before proceeding to 
analyze and fit the experimental data. Finally, in Section 
|VII| we discuss the significance of the present work and 
provide an outlook for future studies. 



where /3 = l/ksT and Dq is the short time diffusion coef- 
ficient at infinite dilution. The time-ordered exponential 
function in Eq. ^ arises because fit (t) does not commute 
with itself for different times . 

An important approximation underlying Eq.([3| (and 
thus p])) is that solvent induced hydrodynamic interac- 
tions (HI) between the colloidal particles are neglected. 
The diffusion coefficient entering Eq.(|3| is thus a scalar 
quantity and the external flow may be included using 
a prescribed (as opposed to self-consistently calculated) 
shear field j{t). While the omission of HI may be inap- 
propriate at high shear rates, for which hydrodynamically 
induced shear thickening can occur in certain systems, it 
is expected to represent a reasonable approximation for 
slowly sheared states close to the glass transition. Nev- 
ertheless, when attempting to fit experimental data us- 
ing theoretical models based on Eq.([3| it proves necces- 
sary to include an empirical hydrodynamic correction ac- 
counting for the high frequency viscosity. In addition to 
the neglect of HI we make two, potentially more dan- 
gerous, assumptions: (i) ^{t) is taken to be spatially 
translationaUy invariant, which may become question- 
able when considering the flow response of dynamically 
arrested states, (ii) The shear field acts instantaneously. 
While this should be acceptable for certain flow histories 
the general status of this approximation is not clear. 

B. Linear response 



II. FUNDAMENTALS 
A. Microscopic starting points 

The shear stress resulting from a general time- 
dependent shear strain of rate "/(t) is given by a gen- 
eralized Green-Kubo relation |T9l [20] 



a(t) 



dt' j{t')G(t,t') 



(1) 



Equation ([!]) is nonlinear in the shear rate due to the non- 
linear functional dependence of the shear modulus G(t, t') 
on 7(i). Within the microscopic framework developed in 
[TOl 12"0"] the modulus is identified as the correlation func- 
tion of fluctuating stresses 



G^t')= } ^ fV {a xy J} dsnHs) a xyi 



(2) 



where <r xy = — J2i FiVi is a fluctuating stress tensor el- 
ement, formed by a weighted sum of the forces acting 
on the particles for a given configuration, T is the tem- 
perature, V is the system volume and (•) indicates an 
equilibrium average. The particle dynamics to be con- 
sidered in the present work are generated by the adjoint 
Smoluchowski operator [3T] 



fit (t) = 2 Do [ S« + PFi }-d l + D 7(*)W 



8 
dx 



(3) 



Following standard convention, we consider an exter- 
nally applied shear strain of the form 



7(<) = 7osin(a;t). 



(4) 



The time translational invariance of the shear field (4) 
gives rise to an explicit dependence of the modulus (2) 
upon two time arguments. 

For small deformation amplitudes (70 <C 1) the strain 
dependence of the shear modulus may be neglected, such 
that Eq.Q provides a linear relationship between j(t) 
and cr(t). This leads to the approximation 



G{t,t') = G e(l {t-t l ) : 



(5) 



where G eq (t) denotes the time translationaUy invariant 
equilibrium shear modulus. Substitution of Q and ^ 
into Eq.Q and employing trigonometric addition formu- 
las leads directly to the familiar linear response result 

a(t) = j G'(u}) sin(cjt) + 7 G"(u;) cos(arf), (6) 

where G'(oj) and G"(ui) are the storage and loss moduli, 
respectively, defined by 



G'(ui) = ui \ dt! sin(tjf') G cq (t') 



G"(uj) = u dt'cos(ut')G ecL (t'). 



(7) 
(8) 
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Furthermore, Eq.([6]) can be rewritten as 

<7(t)=7b|G(w)|sin(wt + 5(a;)), (9) 

where the complex modulus is given by G = G + iG" 
and the phase shift by S = arctan(G"/G). If G"{u) = 
the response is purely elastic, in phase with j(t) (5 = 0). 
In the case G'(ui) — dissipation dominates and the 
response is in phase with j(t) (5 = 90°). 

C. Nonlinear response 

It should be clear at this stage that the familiar linear 
response form ^ is a direct consequence of the convolu- 
tion integral which results from inserting the time trans- 
lationally invariant equilibrium function ^ into Eq.([2]). 
For finite strain amplitudes, the dependence of the modu- 
lus upon two time arguments prevents the simple trigono- 
metric manipulations leading to Eq.([6]). Nevertheless, 
the non-sinusoidal stress response, cr(t), is periodic with 
period 2tt/uj, and may therefore be expressed as a Fourier 
series 



<K*) = 7o ^2 G' n {uj) sm(nujt) + 7o X! G "M cos(nujt), 



n=0 



(10) 



where G' n and G" are frequency dependent Fourier coef- 
ficients given by [H] 



ir/u 



dt a(t) sin(no;t) 



— 7r/ Id 
7r/u) 



G"(w) = - / dt a(t)cos(nut). 



(11) 
(12) 



In the limit 70 — > the coefficients G[ and G'{ reduce to 
the familiar linear response moduli. It should be noted 



that we retain the n — term in the second sum of ( 10 ) 



in order to leave open the possibility of a stress offset. 

Employing manipulations analogous to those leading 
from Q to ([9| the Fourier series ( 10 1 may be expressed 
in the form 



c(t) = 7o sin(nwt + 5 n (cj)), 

n=l 



(13) 



where the amplitude is given by I n — \G' n + iG''| and the 
phase shifts by 5 n (oj) = arctan(G"/G^). In analyzing our 
theoretical, experimental and simulation results we will 
focus on the behaviour of both the generalized moduli G' n 
and G" and the amplitude and phase shift, /„ and 8 n , of 
the fundamental (n = 1) and higher harmonics (71 > 1) 
as a function of the control parameters. 

Following a period of transient response after initia- 
tion of the strain field (switching on the rheometer) the 
system enters a stationary state, demonstrating a peri- 
odic stress response. Although, to some extent, an is- 
sue of semantics, it is important that the 'stationary' 



state presently under consideration be distinguished from 
'steady' states, of the kind achievable by application of 
a time-independent shear flow. The stationary state is 
simply a well characterized and periodic transient and is 
thus influenced by additional physical mechanisms (e.g. 
thixotropy) which are irrelevant for steady states. In a 
physical system the stationary response must be inde- 
pendent of the direction of shear, leading to a stress a(t) 
symmetric in j(t). The mirror symmetry of the consti- 
tutive equation has the consequence that only odd terms 
contribute to the Fourier series (13 1. The appearance of 



even harmonics in the analysis of experimental data is 
often an indication of boundary effects, such as wall slip, 
or other inhomogeneities of the flow }10j . 

Important physical interpretation may be given to the 
coefficient G'[ by considering the energy dissipated per 
unit volume of material per oscillation cycle 



7r/a; 

£',/ = / dt<r{t)j(t). 

-7T /U) 



(14) 



Substitution of the strain field Q and the Fourier series 
(fl0| into Eq.JII} leads to 



E d = ^G'l{w), 



(15) 



(see also [H]). Thus, for a sinusoidal strain of the form 
Q, energy is dissipated only at the input frequency. The 
coefficient G" therefore has the same interpretation in the 
non-linear regime as in the linear regime: it determines 
the dissipation of energy over an oscillation cycle. The 
remaining coefficients in the series, G' n and G" >1; thus 
collectively describe the reversible storage and recovery 
of elastic energy. 



D. Lissajous plots 

A standard way to graphically represent the relation- 
ship between j(t) and a(t) is via the Lissajous repre- 
sentation, in which trajectories are shown in the 7*,er* 
plane, where 7* = 7/7 max and a* = cr/cr m ax are the 
strain and stress, normalized by their maximum values 
[3 7) . In this representation, a general linear viscoelastic 
response is characterized by an ellipse, symmetric about 
the line 7* — a*, point symmetric with respect to the 
origin plus two mirror planes. The two limiting cases of 
a purely elastic and a purely dissipative response are thus 
characterized by a line and a circle, respectively. In the 
nonlinear regime considerable deviations from ellipticity 
are observed. The specific character of these deviations 
can indicate whether a material is, for example, strain 
hardening or strain softening (an increase/decrease of G' 
with strain amplitude), and thus provides a useful, al- 
beit qualitative, 'rheological fingerprint' of a given ma- 
terial [6j|8j. For a general nonlinear response, the area 
enclosed within the closed loop trajectory of a Lissajous 
figure is directly related to the dissipated energy via the 
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integral in Eq. (14 1. This lends an appealing physical in- 



terpretation to the Lissajous representation and provides 
a direct visual impression of the dissipative character of 
the response. 



III. THEORETICAL APPROACH 
A. Schematic model 

As noted in the introduction, the approximate micro- 
scopic constitutive theory developed in [THl HO] enables 
first-principles prediction of the rheological behaviour of 
dense colloidal dispersions. However, the simultaneous 
occurrence of spatial anisotropy and non-time transla- 
tional invariance hinders numerical solution of the equa- 
tions when addressing concrete problems. The schematic 
model presented in 24J provides a simplified set of equa- 
tions which, it is hoped, capture the essential physics 
contained within the full equations, while remaining 
tractable for numerical implementation. 

Within the schematic reduction, the modulus is ex- 
pressed in terms of a single-mode transient density cor- 
relator 



G{t,t') = v a <S> 2 (t,t') 



(16) 



where v„ is a parameter measuring the strength of stress 
fluctuations. The approximation underlying ( 16 ) is that 



stress fluctuations relax as a result of relaxations in the 
density (viz. structural relaxation) . The assumption that 
v a is independent of strain is a simplifying assumption 
which could be relaxed if neccessary. The microscopic 
theory of [2D] predicts both the temporal and wavevec- 
tor dependence of the transient density correlater under 
applied flow. The schematic, single mode, density cor- 
relator (normalized to $(£,£) = 1) represents, in some 
non-specific sense, a 'typical' correlator of the micro- 
scopic theory. It is obtained from solution of a nonlinear 
intcgro-diffcrcntial equation 

<f>(t,t')+r(<f>(t,t')+ [ dsm{t,s,t')<Z>(s,t')) = 0. 



(17) 

The single decay rate T sets the time-scale and would, 
within a microscopically based theory, depend upon both 
structural and hydrodynamic correlations. The overdots 
in Eq.(17) imply differentiation with respect to the first 
time argument. The memory function m(t, s,t') appears 



in Eq.(17) as a generalized friction kernel, which can be 
formally identified as the correlation function of certain 
stress fluctuations. Making the assumption that these 
stress fluctuations may be expressed in terms of density 
fluctuations (both become slow close to the glass tran- 
sition) leads to a tractable expression for m(t, s,t') as a 
quadratic functional of the transient density correlator 
and, thus, a closed theory. A somewhat surprising con- 
sequence of the formal calculations presented in [TO] [20] 



is that the memory function possesses three time argu- 
ments. The presence of a third time argument, which 
would have been difficult to anticipate on the basis of 
quiescent MCT intuition, has important consequences for 
rapidly varying flows (e.g. step strain [11]). Within the 
schematic model the memory function is given by 



m(t, s, t') = h(t, t') h(t, s) [vx$(t, s) + v 2 ^ 2 {t, s) ' 



(18) 



Following conventional MCT practice the parameters V\ 
and v 2 represent, in an unspecified way, the role of the 
static structure factor in the microscopic theory and are 
chosen as v 2 = 2 and v x = 2(%/2 - 1) + e/(\/2 - 1). The 
separation parameter e is a crucial parameter within our 
approach and encodes the thermodynamic statepoint of 
the system by measuring the distance from the glass tran- 
sition. Negative values of e correspond to fluid states and 
positive values to glass states. Setting h equal to unity 
in Eq.(18l recovers the well known F\ 2 model, originally 



introduced by Gotze |26l 03] [H] . The linear term in <f> 
which appears in Eq.(18l is absent from the microscopic 



mode-coupling expression, but turns out to be necessary 
for a faithful reproduction of its asymptotic properties 
within a single-mode theory. Under simple shear flow, 
the /i-functions in the memory kernel ( 18 1 serve to accel- 



erate the loss of memory caused by the affine advection 
of density fluctuations. The assumption that the same 
function h may be used to incorporate both the strain 
accumulated between t and t' as well as that between t 
and s is an approximation, made to keep the theory as 
simple as possible. Taking account of the required invari- 
ance with respect to flow direction suggests the simple 
ansatz 



h(t,t') 



lc+l 2 ' 



(19) 



where 7 = ^(t, t') — f., ds 7(s) and the parameter 7 C sets 
the scale of strain. 

Eq.Q and fl6])-([l9]) provide a closed constitutive the- 
ory which depends upon three adjustable parameters 
(f CT ,r,7 c ) and two control parameters (£,7) represent- 
ing the coupling strength and applied shear rate. As the 
schematic model under discussion is implicitly based on 
the Smoluchowski dynamics described by Eq.([3|, the in- 
fluence of HI is neglected. While this is not important for 
capturing correctly the qualitative features of the rheo- 
logical response, quantitative comparison requires a sim- 
ple hydrodynamic correction at high frequencies. The 
simplest approximation, which we will employ in the 
present work, is to empirically add an extra term to the 
shear modulus 



G{t,t')^G{t,t') + r} 00 8{t-t'). 



(20) 



The high frequency viscosity, rj^ , is thus introduced into 
the model, describing the viscous processes which occur 
on timescales much shorter than the structural relaxation 
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time. The correction (20) has the consequence that the which, together with the equation of motion 



stress acquires an extra term, linear in 7, and the Fourier 
coefficient G'{ is shifted by a term linear in u). 



B. Strain- rate frequency superposition 

An alternative mode-coupling-type approach, describ- 
ing the collective density fluctuations of dense colloidal 
fluids under shear, is provided by the work of Miyazaki 
et al. [3"5H4Ti] . By considering time-dependent fluctua- 
tions about the steady state a closed (scalar) constitutive 
equation has been derived and applied to colloidal disper- 
sions in two-dimensions under steady shear [381139] and in 
three dimensions (subject to additional isotropic approx- 
imations) under large amplitude oscillatory shear [40] . 
Given the very different nature of the approximations un- 
derlying the present MCT-based theory |2"0"I [2"2"l |2"5] 
and that of [3"5H3U] (fluctuating hydrodynamics vs. pro- 
jection operator methods) it is interesting that the final 
expressions (e.g. the memory function vertices entering 
the equation of motion for the transient correlator) are 
rather similar, at least for the special case of steady shear. 
For the case of large amplitude oscillatory shear, however, 
the theory presented in [40 differs clearly and fundamen- 
tally from the microscopic approaches to time-dependent 
shear developed in [T3 [2D] and, consequently, from the 
schematic model of [24] to be employed in the present 
work. The theoretical developments of Miyazaki et al. 
[iO] motivated the authors to propose the principle of 
'strain-rate frequency superposition' as a probe of struc- 
tural relaxation in soft materials [45] . 

The essence of the Miyazaki et al. approach can be 
captured by a simple schematic model, which we will 
elaborate upon below. In 40J the authors took the 
theory which they had developed for steady shear flow 
[551 I3T?] and replaced the steady shear rate 7 appear- 
ing in the equation of motion for the correlator, by the 
time-dependent shear rate j(t) — ^quj cos(ujt) , describing 
oscillatory flow. This rather ad hoc treatment gives rise 
to equations with a mathematical structure appropriate 
for steady flows and ignores the more realistic, although 
more complicated, history dependence of theories devel- 
oped to treat non-steady flows specifically [13 [23 . On 
the basis of the results obtained for the strain amplitude 
dependence of the storage and loss moduli (notated as 
Gi , G'{ in the present work) it was argued that the time- 
dependence of the strain-rate field 7(i) = 7 wcos(cj£) is 
not essential for understanding the viscoelastic response, 
and that it is sufficient to consider the strain-rate am- 
plitude 70W alone. The relevant timescale is thus iden- 
tified as (7oa;)~ 1 , rather than us^ 1 . Within the context 
of schematic mode-coupling equations, this assumption 
may be expressed by the following memory function 



m(t) 



[v 1 ^{t)+v 2 ^ 2 {t)' 
I + ( 7 (M) 2 



$(t) + r( $(*) + J dt' m(t - t')${t') ) =0, 



the shear modulus 



G(M') = v a <S> 2 {t-t'), 



(22) 



(23) 



and Eq. pT) p rovides a closed theory for <r(t). In fact, 
Eqs.(2l|-|23) are identical to the F? 2 model [23127], with 
a steady shear-rate 7 = 70 u). An important consequence 
of assuming the dominance of the timescale (70 w) _1 is 
that all states, even those which would be glasses in the 
absence of flow, become fiuidized by an applied oscilla- 
tory shear field, regardless of the amplitude 70. Whether 
or not a vanishingly small value of 70 is really sufficient 
to restore ergodicity to dynamically arrested states is un- 
clear and presents a fundamental question, with impor- 
tant implications for the existence of a linear response 
regime. 

Despite capturing approximately the amplitude depen- 
dence of G' 1 ,G'{, describing the response at the funda- 
mental frequency, higher harmonics are ignored in the 
approach of [30]- The absence of higher harmonic con- 
tributions within the theory of Miyazaki et al. can be 
traced back to the assumption that the time-dependence 
of 7(i) is irrelevant and that this can be represented by 
the constant 70^. Within the present context this has 
the consequence that the memory function (21 ) and cor- 



(21) 



relator, given by solution of (22 1, are constrained to be 
time-translationally invariant (viz. depend on a single 
correlation time only). While this assumption is clearly 
at odds with the underlying variations in the strain field 
(for which a dependence on the 'waiting time' is to be ex- 
pected), it nevertheless serves to capture first-order cor- 
rections to linear response theory, while remaining rela- 
tively easy to implement numerically. 

The theory developed in [40] is quasi-linear, in the 
sense that a(t) remains a simple sinusoid, but with an 
amplitude and phase shift which depend non-linearly on 
70. Attempts to justify the neglect of higher harmonics 
have been based on the fact that the ratio of the third 
harmonic amplitude to that of the fundamental remains 
smaller than approximately 20%, for a wide range of sys- 
tems [33 S3- However, in order to draw a fair conclu- 
sion, it is important to consider the sum J2 n >i In(w)/Ii, 
rather than I3/I1 alone, when assessing the physical rel- 
evance of higher harmonic contributions. Various exper- 
imental studies on colloidal dispersions (see e.g. [33]) 
show clearly that the higher harmonics can collectively 
account for up to half of the total signal, which is not 
a small effect. This observation serves to emphasize the 
importance of truly nonlinear theories, which confront 
directly the non-time translational invariance of the cor- 
relation functions, thus going beyond the convolution ap- 
proximation to Eq.Q. 
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IV. COMPUTER SIMULATION 

To provide a point of reference for our theoretical cal- 
culations we have performed two-dimensional simulations 
on a system hard discs undergoing Brownian motion in 
an external shear field. The simulations are designed to 
solve approximately the many-body problem of a system 
of interacting Brownian particles under shear flow. The 
same Smoluchowski dynamics [3T] underlies the micro- 
scopic mode coupling theories of [TH] and [2U] which form 
the basis of the schematic model employed in the present 
work [21] . We choose to simulate a two-dimensional sys- 
tem for two reasons: (i) The computational resources 
required are significantly reduced with respect to sim- 
ulation of three-dimensional systems and thus enables 
improved statistics to be obtained, (ii) Recent micro- 
scopic studies of the quiescent mode-coupling theory in 
two-dimensions have revealed behaviour broadly similar 
to that found in three dimensional calculations [46]. We 
thus expect the reduced dimensionality of our simulation 
system to be of little consequence for qualitative compar- 
ison with the present theory and experimental data. 

The basic concept of the algorithm has been described 
in detail in three dimensions in [55 and its adaptation 
to two dimensions can be found in |56) . We consider 
a binary mixture of hard discs with the diameters of 
D s = 1.0 and Df, = 1.4 with equal particle number con- 
centrations and a total amount of N = N s + Nb = 1000 
hard discs in a two-dimensional simulation box of volume 
V with periodic Lees Edwards boundary conditions. The 
total two-dimensional volume fraction is then given by 
0tot = ^pK-Dft + D 2 S ). We employ this system in order 
to suppress crystallisation effects. The mass m of the 
particles and fc^T is set equal to unity. We choose our 
coordinate axes such that flow is in the x-direction and 
the shear gradient is in the y-direction. The Brownian 
timestep was chosen to St — 0.01 as in [55] . This results 
in a short time diffusion constant of Do = 0.05. To im- 
plement a time-dependent, oscillatory shear rate, at each 
Brownian timestep the shear rate is set to its new value 

7(tb) = 7o wcos(W B ), (24) 

and all particle velocities are freshly drawn from the 
Gaussian distribution with (v 2 ) = 2 and (v) = 
7( T s)2/( T s)- Between two Brownian timesteps the shear 
rate is kept constant. The strain "f(t) can therefore be 
obtained using 

7(*)= E 7(*i0ft, ( 25 ) 
r B e[o,t] 

which leads to 7(t) = L dt j(t) in the limit of St — > 0. At 
every Brownian timestep the part 7(tb)j/(tb) guarantees 
a linear velocity profile as a linear shear flow is imposed 
on every particle, depending on its y-position. For all 
simulations the frequency was set to w = 0.001 which 
leads to Pe w = ujD 2 /AD q = 0.05. 



The average quantity of interest in the present work 
is the time dependent potential part of the shear stress 
<Txy(t) = y (Y,( F ( t )ij)x{r(t) l j) y ) , with the relative force 
components of particle i and j (F(t)ij) x and the particles 
relative distance component (r(t)ij) y for a given time t. 
As we consider hard particles the forces must be calcu- 
lated from the collision events. By observing the colli- 
sions within a certain time window Ar c = [tf. , tk + Ar c ] 
for a given time tk, forces may be extracted using the 
change of momentum which occurs during the observa- 
tion time. This leads to the evaluation algorithm for the 
stress at time tk 

\ c t c e[t k ,t k +AT c ] I s 

where summation is over all collisions after time t c 
within the time window Ar c . The procedure effectively 
sums the momentum changes (Avij) x in x— direction 
multiplied by the relative distance of the particles (rij)y 
in y— direction. The brackets (...) s denote the different 
simulation runs. 

At a total volume fraction of <f>tot = 0.81 which is slightly 
above the glass transition for this system (estimated to be 
at ip = 0.79 on the basis of simulated flow curves [SS]) we 
prepared 4000 independent sets for each amplitude 70 € 
{0.001, 0.003, 0.009, 0.01, 0.09, 0.1, 0.2, 1.0, 10.0, 100.0} 
As the system starts from a non stationary state it 
is necessary to wait for the system to reach it's long 
time asymptote (which we found to be the case after 
undergoing two full oscillations) before meaningful 
averages can be taken. 



V. EXPERIMENT 

A. Characterization of the latex particles 

The polydisperse latex particles consist of a solid core 
of poly(styrene) (PS) onto which a thermosentitive net- 
work of crosslinked poly(N-isopropylacrylamide) (PNI- 
PAM) is affixed [35]. The degree of crosslinking of the 
shell due to the crosslinker N,N'-methylenebisacrylamide 
is 2.5 Mol %. As exactly the same particles were used for 
this work as in |32| . the latices have a temperature de- 
pendent size (hydrodynamic radius in nm Ru — -0.7796 
•T +102.4096 with T the temperature in °C below 25°C) 
and a polydispersity of 17% 32 . All experiments were 
done in an aqueous solution of 0.05M KC1 to screen resid- 
ual charges which emerge from the synthesis of the parti- 
cles. The solid content of the suspension was determined 
by comparing the weight before and after drying and was 
found to be 8.35 wt%. Neither the MCR 301 measure- 
ments nor the FT-rheology measurements at 15.1°C lead 
to a significant change of the solid content (+0.02 wt%). 
However, the FT-rheology measurements at the remain- 
ing two temperatures (18.4°C and 20°C) had a slightly 
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different solid content (9.02 wt%) due to the physical 
relocation of the rheometer to another laboratory and 
some additional solvent evaporation. The effective vol- 
ume fraction (j) e ff was calculated by using the correlation 
of mass concentration c, hydrodynamic radius and effec- 
tive volume fraction found in the inset of fig. 6 in |32j . 
which is given by c • Ffi H = 9.67 • 10~ 17 <? • (f> e s for the 
different temperatures. For the temperature of 15.1°C a 
volume fraction <f) c g of 0.65, for 18.4° C a <t> c s of 0.60 and 
for 20.0°C a o ff of 0.57 was found. In previous work [32] 
the glass transition for this system was found to be at 
4>cS = 0.64. Given a polydispersity of 17% the theory of 
Schaxtl and Sillescu [33] predicts random close packing 
at 4> e ff = 0.68. 



Rheometer 

r.^r,^, ^ u , ^ (Rheometric Scientific 

PC RSI-Orchestrator 




FIG. 1: (Experiment) The experimental set-up of FT- 
rheology with two independent controlling computers, one for 
normal rheology and one for FT-rheology [5]. 



B. Rheological Experiments 

The rheological experiments cover the range from the 
linear to the strongly nonlinear regime. The rheological 
measurements were performed with two instruments a 
MCR 301 from Anton Paar with a Cone Plate geom- 
etry (diameter: 50 mm, cone angle: 0.991°) and the 
FT-rheological measurements with an ARES rheometer 
(Rheometrics Scientific) with a Cone Plate geometry (di- 
ameter: 50 mm, cone angle: 0.04 rad). All measurements 
were performed after a preshear of 100s _1 lasting 200s 
and a waiting time of 10s. For the rheological measure- 
ments, a solvent trap (for the ARES instrument equipped 
with a sponge drawn with water) is used to prevent evap- 
oration. A thin paraffin-layer covers additionally the so- 
lution to prevent an exchange with the atmosphere above. 
With the MCR 301 the measurements at 15, 18 and 20°C 
were done. The ffowcurves were measured from 5-10 -4 
to 1000s -1 with a time ramp of 2500s to 20s and then 
a waiting time of 10s followed by a ffowcurve measured 
from 1000 to 5-10~ 4 s _1 with a time ramp of 20 to 2500s. 
For 15° C ffowcurves in the shear rate range of l-lO^s -1 
to 1000s -1 with a timeramp of 10000s to 20s were per- 
fomed. The frequency tests were performed at a defor- 
mation of 1% starting from 10Hz to 0.001Hz with a time 
ramp of 20s to 1000s. All deformation tests were per- 
formed with a measurement time for each point of 100s. 
For the FT-rheology oscillatory time sweep measure- 
ments were performed at the ARES instrument at fre- 
quencies of 1, 0.1 and 0.01Hz at different deformation 
amplitudes. The FT-signal were always recorded after 
some oscillations, so that the suspension reached the os- 
cillatory stationary state. The first measurements were 
done at 15.1°C and the frequency tests of the ARES and 
the MCR 301 coincided. For the two higher tempera- 
tures the measurements had to be performed for the same 
ARES instrument at a different room. Here an evapora- 
tion effect was observed, so that the temperature had to 
be adjusted to archieve the same behaviour in the fre- 
quency test. Instead of 18°C a temperature of 18.4°C 
and instead of 20°C a temperature of 20.9°C was used. 
Typically for the nonlinear FT-rheology measurements 



with the time sweep tests 40 oscillations for 1Hz excita- 
tion were applied, whereas 10 oscillations for 0.1Hz and 
9 oscillations for 0.01Hz. To obtain a FT-spectrum from 
the raw time data we performed a discrete, complex, half- 
sided fast Fourier transformation [5] I12[ 154] . 

For further information of the setup, the measuring 
principal and FT-analysis we would like to refer to |51[ 

[S2] 



C. Fourier Transform Rheology Apparatus 

The FT-rheological setup consists of a Rheomet- 
rics Scientific advanced rheometer expansion systems 
(ARES) and a computer which either controls the 
rheometer via a serial cable as well as detects the strain 
and torque force outputs via BNC cables (see Fig[T]). 
The ARES rheometer is a strain-controlled rheometer 
equipped with a dual range force rebalance transducer 
(100 FRT) capable of measuring torques ranging from 
0.004mNm to lOmNm, specified by the manufacturer. It 
has a high resolution motor, applying frequencies from 
10 — 5rad/s to 500rad/s and deformation amplitudes 
ranging from 0.005 to 500mrad. A water bath adjusts 
the temperature range from 5°C to 95° C. The analog 
raw data of the measurements are digitized with a 16-bit 
ADC. This ADC card has a maximum sampling rate of 50 
kHz per channel. Due to the high sampling rate the time 
between consecutive data points is very small compared 
to the timescale of rheological experiments. The loss of 
information by sampling the torque transducer data is 
negligible [48 . To achieve best results with respect to 
the signal-to-noise ratio, oversampling is applied. The 
ADC-card acquires the time data at the highest possible 
sampling rate and then preaverages them on the fly to 
reduce random noise. With this method the noise is re- 
duced by a factor of 3 to 5 which could only be achieved 
by averaging multiple measurements |49j . Within the 
set-up a 16-bit ADC card is implemented, which is able 
to discriminate steps. The quantification resolution of 
the ADC card limits the ratio. It determines the mini- 
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mum detectable intensity of weak signals by its ability to 
discriminate the intensity of the signal. The higher the 
bit number, the smaller the detectable intensity variation 
[50) . After acquisition and digitization of the time data, 
they are handled with Matlab software [52J • 



VI. RESULTS 
A. Theoretical predictions 

1. Flow curves 

For given values of the parameters (« g ,7 c ,r,e) the 
schematic theory defined by ([I]) and Eqs.([ll3|)-([T9"| en- 
ables prediction of the flow curve expressing the steady 
shear stress a as a function of shear-rate j. Fig|2] shows 
a set of typical flow curves generated by the schematic 
MCT model for three fluid states (e < 0), the critical 
state (e = 0), and three states in the glass (e > 0). The 
parameters employed for the theoretical calculations pre- 
sented in Fig(2] (as well as for Figs 3][S I are v a = 100, 
7 C = 0.15, T = 100 and e = 10~ 3 . Experience with fit- 
ting the experimental data, to be considered in section 
|VI C[ shows that these choices represent sensible physi- 
cal values for the model parameters. In the fluid, there 
exists a linear (Newtonian) regime for small shear-rates 
(Pe = 7/r <C 1), for which the standard F 12 model re- 
sult for the shear viscosity holds, a = 7 / dt&l q (t) = 
777. Increasing the separation parameter to less negative 
values (corresponding to, e.g. an increase in the volume 
fraction) gives rise to an increase in 77, reflecting the slow- 
ing of the structural relaxation time r, which dominates 
all transport properties within our MCT approach. For 
7T > 1 the effect of shear starts to dominate the struc- 
tural relaxation and the stress increases sub-linearly as a 
function of shear-rate, corresponding to shear thinning of 
the viscosity 77(7) — a^/j. At high shear-rates jr ^> 1 
the present model yields a = j/T and needs to be supple- 
mented by corrections which account for the high shear 
limiting viscosity (and which, in the absence of HI, are- 
determined by the solvent contribution 7/00). 

As e —> 0~ the regime of linear response shifts to in- 
creasingly lower values of the shear-rate and disappears 
entirely at the (ideal) glass transition, e = 0. For states 
in the glass there exists a finite stress in the limit of van- 
ishing shear rate, identified as the dynamical yield stress 
(linXy^o 17(7) = (J y ). Within idealized MCT based treat- 



ments the dynamical yield stress emerges discontinuously 
as e is varied accross the glass transition (shown in the 
inset of Fig{2]). It should be mentioned that the flow 
curves shown in Fig[2j differ quantitatively from those of 
the extensively studied F± 2 model [27], due to the inclu- 
sion of an additional prefactor h(t, t') in the expression 
for the memory function (18). Nevertheless, the qualita- 
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FIG. 2: (Theory) The flow curves for fluid states e = 
— 10~ 3 , —5 x 10~ 4 , — 10~ 4 (green), at the critical point e = 
(red) and for glassy states e = 10~ 4 ,5 x 10" 4 , 10~ 3 (blue). 
In the glass (e > 0) there exists a finite stress in the limit 
of vanishing shear rate, identified with the dynamical yield 
stress (lim-y-x} = a v ). The cross indicates the yield stress 
value used in Figs[4]and[5] Inset (a) shows the discontinuous 
emergence of a dynamical yield stress as a function of e. Inset 
(b) shows the viscosity r\ = a/7. Calculations were performed 
with ?7oc = 0. 



2. Linear response moduli 

The linear storage and loss moduli, given by Eqs.Q 
and ^ , respectively, are shown in Fig{3] as a function 
of Pe w = uj/T for two fluid states (e < 0) and two 
glassy states (e > 0). In the fluid, the finite value of 
the structural relaxation timescale r is reflected in the 
maximum of G" and the crossing of G' and G" at low 
frequency. The fact that G' remains notably larger than 
G" at high frequencies is simply a result of neglecting 
the high frequency limiting viscosity in presenting 
our theoretical predictions. Setting 7700 = in presenting 
the theory highlights the contribution of structural pro- 
cesses to the viscoelasticity. In the glass, G" goes to zero 
at low frequencies (G" ~ to) and G' attains a finite low fre- 
quency value, identifying the transverse elastic constant 
Goo. Within our MCT approach, the elastic constant 
appears discontinously upon crossing the glass transition 
(i.e. lim^-^o G'(lS) jumps from zero for e — > 0~ to a fi- 
nite value for e — > + ) thus demonstrating that the MCT 
indeed describes a transition to an amorphous solid. 



3. Nonlinear stress response 



By numerical solution of the equation of motion (17) 



tive predictions of the theory for the flow curves are in 
full agreement with those of the F^ 2 model. 



we obtain the non-time translational invariant density 
correlator &(t, t') and, via Eqs.Q and ( 16 ), the nonlinear 
stress response. The numerical algorithm requires the 
equation of motion to be discretized over the entire two- 
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FIG. 3: (Theory) The linear response moduli G[ (full lines) 
and G" (broken lines) as a function of Pe^ for two statepoints 
in the fluid e = -0.001 (red), -0.0005 (green) and two in the 
glass e = 0.0005 (blue), 0.001 (black). For fluid states the 
finite value of the structural relaxation timescale is reflected 
in the maximum in G" and the consequent crossing of G' 
and G" at low frequency. The results presented here omit 
solvent hydrodynamics which may become relevant at Pe^ ~ 
1 (calculations were performed with r/^, — 0). These will 
be included in a simple approximate fashion when fitting the 
experimental data. The two arrows indicate the values Pe w = 
0.001 and 0.025 used to generate Figs|4]and[5j respectively. 



dimensional (t, i')-plane. While considerations of both 
causality (i > t 1 ) and the periodicity of the correlator 
with respect to a translation in time ($(< + t , t' + t ) = 
&(t,t'), where to = tt/uj) enable certain simplifications to 
be made, calculation of the correlator over many decades 
in time remains a computationally demanding task. 

Typical examples of the response a(t) for a glassy state 
(e = 0.001) are shown in Figji] for various values of the 
strain amplitude 70 at a fixed frequency. Alongside each 
of the time series we show the corresponding Lissajous 
curves indicating the extent of the dissipation (via the 
area enclosed, see Eq.(14)) and the deviations from non- 



linearity (discernablc from the non-ellipticity of the loop) . 
The value of e employed to generate this figure generates 
a dynamic yield stress a y = 0.2763, which is indicated in 
each panel of Figj4] by a broken red line. This can also 
be read-off from the appropriate flow curve in Fig(2] 

For small amplitudes 70 < 0.01 the system is almost 
linear and responds in a predominately elastic fashion at 
the considered frequency. As 70 is increased, clear devia- 
tions from a sinusoidal response are apparent and higher 
harmonics start to contribute to the signal. In this non- 
linear regime the stress response exhibits a characteristi- 
cally flattened peak, with an asymmetry which increases 
as a function of 70. It is clear from the figure that the 
higher harmonics first become significant when the max- 
imum value of the stress a(t) approaches the dynamic 
yield stress. An idealized yield stress material, subject to 
large amplitude oscillatory shear of vanishing frequency, 
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FIG. 4: (Theory) The stress response of a glassy state to os- 
cillatory strain calculated from our MCT based theory for 
strain amplitudes from 70 = 0.01,0.03,0.05,0.07,0.10,0.15 
and 0.20. The associated Lissajous figures illustrate the non- 
linear character of the response. The increase in dissipation 
with increasing 70 is reflected by the increasing area enclosed 
by the Lissajous curves. All calculations were performed 
at Pe„ = 0.025 and e = 0.001. The red horizontal bro- 
ken lines indicate the dynamic yield stress obtained from the 
flow curve in Figurej2]for e = 0.001 {a y = 0.2763). The re- 
sponse becomes clearly nonlinear when the maximum of a(t) 
approaches the dynamical yield stress. The blue horizonal 
dotted lines provide an upper bound for the maximum of the 
time-dependent stress and are taken from the corresponding 
flowcurve in Figj2] 



would be expected to show a steady increase of the stress 
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cot 

FIG. 5: (Theory) As in Figure|3] (e = 0.001) but for a lower 
frequency Pe„ =0.001. At this value of Pe^ the system is 
almost perfectly elastic in the linear regime (G' 3> G" for 
70 <C 1, see Figurels). As the time-dependent stress exceeds 
the dynamical yield stress the signal becomes clipped. At this 
frequency cr max lies very close to a y and has thus been omitted 
for clarity. 



up to the yield point, beyond which the system begins to 
flow, maintaining a(t) — o y until the reversal of strain en- 
ables relaxation back to zero. If this were the case, then 
cr(t) would be represented by a 'clipped' signal, symmet- 
ric during loading and unloading of the sample. 

The results presented in Fig[2J however, have been gen- 
erated for a low, but not vanishingly small, frequency. At 
finite frequency the maximum stress attainable in a sys- 



2p 

Nonlinear Response 

Pe =0.025 




L °gio^o 

FIG. 6: (Theory) The theoretical G'i and G" as a function of 
strain for e = 0.001 at frequencies Pe„ = 0.025 and Pe„ = 
0.001. Points are the numerically calculated data points, lines 
are a guide for the eye. For large values of the strain the 
numerical data are well fitted by the power laws G" ~ f^" 
and G' ~ To" 2 " with v — 0.65 (indicated by dotted lines). 

tern under steady shear is given by <7 max = a(jQLu), where 
or (7) is the stress on the steady shear flow curve. This 
maximal stress under flow, <7 max , is indicated by a blue 
dotted line in Fig[4j for the four largest strain ampli- 
tudes considered. In each case, once the stress exceeds 
the yield point the curve flattens, exhibiting a maximum 
which remains bounded from above by cr max . In the low 
frequency limit the lower and upper bounds to the peak 
value of a(t) become equal, a max = a y , such that the sig- 
nal becomes clipped at the yield point. In order to test 
this hypothesis further, we show in Fig(5] stress responses 
generated using the same parameter set as employed in 
Fig|iJ but for a frequency one order of magnitude lower. 
At this reduced frequency, the clipping of a(t) at yield 
is quite clear, although the the peak stress still slightly 
exceeds a y , due to the fact that the values of the strain- 
rate amplitude are not sufficiently small that er max has 
saturated to <j y . 

Two additional comments are in order regarding the 
results shown in Figs|4]and[5] Firstly, for low frequencies 
(e.g. the data shown in Fig{5]) the full time-dependent 
stress signal can be rather faithfully reproduced by the 
simple approximation 

I 7o[G'i(w) sin(wi) + G"(w) cos(wi) ] a < a y 
a{t) = \0 *>a y 

(26) 

where G[ and are the lowest order coefficients in the 
Fourier series (not. The naive picture sketched above is 
thus not completely correct. In order to describe cor- 
rectly the sub-yield response, a(t) < a y , it is neccessary 
to incorporate the 70 dependence of the lowest order co- 
efficients; linear response is insufficient. It is also note- 
worthy that it is the dynamic, yield stress, which plays 
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FIG. 7: (Theory) The theoretical intensities of the third, fifth, 
seventh and ninth harmonic which contribute to the nonlin- 
ear stress response shown in Figure[3] At strains approach- 
ing unity the collective contribution of the higher harmon- 
ics accounts for approximately 50% of the total a(t) signal, 
(e = 10~ 3 and Pe„ = 0.025). 



the crucial role in determining the time-dependent cr(t). 
While the importance of dynamic yield in determining 
the oscillatory response is clear within the present ap- 
proach, it remains to be seen whether this is a constraint 
introduced by employing a prescribed strain or, more sig- 
nificantly, an indication that the dynamic and static yield 
stresses are identical within our approximate theory. The 
simple approximation ( 26 ) contains higher harmonic con- 



tributions as a result of the yield stress clipping criterion. 



4- Fourier analysis 

In order to provide a more systematic analysis of the 
time signal cr(t), we now consider its decomposition into 
Fourier modes and investigate the behaviour of the coeffi- 
cients entering the series ( 10 1 and ( 13 ) as a function of 70 . 



We address first the strain amplitude dependence of G' ± 
and G", thus mimicing the ubiquitous 'strain sweep' ex- 
periments generally used to assess the nonlinear response 
of a given material. In Fig(6]we show typical results for 
the lowest order coefficients as a function of strain, for 
two different values of the excitation frequency. 

For small values of 70 linear response is recovered and 
the values of G[ and G'[ may be read-off from the data 
shown in Fig|3] The linear response regime persists up 
to around 70 = 0.01, beyond which G'[ begins to in- 
crease gradually, reaching a maximum value at around 
70 = 0.1. As noted in subsection II C the coefficient G' x 
is proportional to the amount of energy dissipated per 
oscillation cycle. The increase in dissipation observed 
over the range 70 = 0.01 — > 0.1 is probably connected to 
the increasing disruption of the microscopic 'cage' struc- 
ture of dense glassy systems, induced by the externally 
applied strain field. However, such microscopic interpre- 
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FIG. 8: (Theory) The theoretical values for the first, third 
and fifth phase shift which contribute to the nonlinear stress 
response shown in Figure[3] It should be noted that the 
phases are only physically meaningful, i.e. contribute sig- 
nificantly to the total a(t) signal, when the corresponding 
intensities shown in Figurej7] are non-negligable. (e = 10~ 3 
and Pe^ = 0.025). 



tations remain purely speculative within the present con- 
text of schematic model calculations, for which there is 
no explicit spatial resolution of correlated density fluctu- 
ations. 

Deeper insight into the microscopic mechanisms under- 
lying the observed macroscopic response would be pro- 
vided by solution of the full equations presented in [20] . 
The increase in G'{ is associated with a decrease in G[ as 
a function of 70. In contrast to G", there exists no simple 
physical interpretation of the coefficient G^ in the non- 
linear regime. For strain amplitudes exceeding 70 = 0.01 
the recoverable elastic energy becomes distributed over 
G'± and the higher harmonics, such that G\ loses its spe- 
cial status as a 'storage modulus'. 

For 70 > 0.1 the dissipation G" becomes larger than 
G[, indicating a crossover from predominately elastic to 
predominately viscous response, and both functions ex- 
hibit a monotonic decay. For values of the strain ampli- 
tude larger than unity a regime of asymptotic decay is 
entered, characterized by a well defined power law depen- 
dence. We find that the numerically generated data are 
well fitted by power laws G'{ ~ r )^ u and G\ ~ Jq U , with 
v = 0.65(±0.02) and v' — 2v. Moreover, calculations 
performed at various frequencies show that the exponent 
values are independent 

of uj and thus seem to represent a universal aspect of 
the asymptotic decay within the schematic model. While 
the numerical findings are suggestive of a universal ex- 
ponent 7v, analytical calulation of its precise value has so 
far proved elusive. The primary difficulty in extracting v 
from the theory is that, even in the asymptotic regime, 
the correlator $(i, t') retains a residual dependence upon 
the waiting time t' which does not yield readily to an- 
alytic treatment. The integral for the stress ([l} thus 
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Nonlinear Response 
Pe =0.025 



FIG. 9: (Simulation) The stress response measured in Brow- 
nian dynamics simulations of a binary hard-disc mixture un- 
der oscillatory flow. A size ratio 1 : 1.4 was used to supress 
crystallization. The considered strain amplitudes range from 
70 = 0.01 to 70 = 10. The left column of figures shows the as- 
sociated Lissajous curves illustrating the nonlinear character 
of the response. The simulations are performed at <j>tot = 0.81 
(slightly beyond the glass transition, according to our simu- 
lation estimates). The Peclet number is Pe^ = 0.05). 



consists of a complicated superposition of correlators for 
different values of t' . 

An important numerical prediction of the schematic 
model is that the exponents dictating the decay of G 
and G" are related, to within numerical accuracy, by a 
factor of 2. The relation v' = 2v has been observed in ex- 
perimental studies of a variety of soft materials (see sub- 
section 
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FIG. 10: (Simulation) The simulation G[ and G" as a func- 
tion of strain for two dimensional volume fraction 0.81 and 
Peclet number Pe^ = 0.05. For large values of 70 the moduli 



obey the power laws G" 



7o 



and G[ ~ j 



relationship v' — 2v is found also in simple nonlinear 
Maxwell models [13 HO] , and is thus not particularly sur- 
prising. Such models inevitably predict the trivial expo- 
nents v = 1 and v' = 2. The numerical data obtained by 
Miyazaki et al. [10] from solving their microscopic MCT 
theory is consistent with the exponent relation v' = 2v, 
but predicts a value of v = 0.9 (obtained by fitting the 
numerical data) , which differs somewhat from the exper- 
imental value v — 0.7 for PMMA colloids presented in 
the same work. Whether the value v = 0.9 from the 
MCT calculations of [40] is influenced by the additional 
isotropic approximations employed remains unclear. 

The coefficients G^ and G'{ discussed above describe 
the response at the fundamental frequency. We now con- 
sider the contribution of the higher harmonic terms to 
the stress signal. Due to the j(t) — — j(t) symmetry of 
o~(t), even coefficients in the Fourier series ( 13 ) are identi- 



VI C for more details on this point) and it would 



therefore be of considerable interest to investigate this 
apparent prediction of our model in more depth. The 



cally zero within the schematic theory (a condition which 
provides a useful check for our numerical algorithms) . In 
FigJT] we show the intensities of the odd harmonics (nor- 
malized by 1\ ) up to n = 9 for a glassy state, obtained by 
applying a discrete Fourier transform to the time series 
a(t). For very small amplitudes 70 < 0.01 the numerical 
solution of the equation of motion for §(t,t') becomes 
unreliable, as the structural relaxation time exceeds the 
range of the numerical 

grid upon which the oscillating function $(t, t') can be 
resolved. Data are thus presented for 70 > 0.01 where 
accurate converged solutions can be obtained. At strain 
amplitude 70 = 0.01, only ^3 contributes significantly to 
the signal (around 3%). As the 70 is increased beyond 
0.01, the increasing influence of 73 is accompanied by the 
appearance of terms ^5, I7 (beyond 70 = 0.03) and Ig 
(beyond 70 = 0.07). Although intuitive, it is not clear 
a priori that the higher harmonics must neccessarily ap- 
pear in sequence n = 3, 5, 7, • • • upon increasing the am- 
plitude. All of the I n >3 exhibit a maximum in the range 
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FIG. 11: (Simulation) The normalized intensities of the third FIG - 12: (Simulation) The phase shifts of the first and third 
and fifth harmonic contributing to the nonlinear stress re- harmonic contributing to the nonlinear stress response shown 
sponse shown in Figure[9] m Fig ure |9] 



0.3 < 70 < 1 and by 70 = 1 contribute approximately 
half of the total signal, Y^ n >i In/h ~ 0.5. The maxi- 
mum and subsequent decay of the higher harmonics has 
also been observed in [S3] . 

Complementary to the higher harmonic intensities are 
the phase shifts S n shown in FigjHJ It should be noted 
that the phases are only physically meaningful for am- 
plitudes at which the corresponding intensity is signifi- 
cant. As 70 is increased towards unity the phases satu- 
rate to the asymptotic values 5± = 7r/2,<53 = 37r/2 and 
85 = 5tt/2. The higher harmonic contributions /„ and 5 n , 
which contribute for 70 > 0.01, describe the distortion of 
a(t) close to the yield stress (c.f. Figsj4]and[5]). 



B. Simulation results 

In Fig|9]we show the stress response measured in our 
Brownian dynamics simulations of a binary hard-disc 
mixture. As the strain amplitude is increased the simu- 
lated stress evolves from a linear to a nonlinear response 
for 70 > 0.03. Consistent with the theoretical results 
shown in Fig|4j the time dependent signal becomes dis- 
torted away from a pure sinusoid when the peak of a(t) 
encounters the dynamic yield stress. Although the Pcclct 
number Pe^ — 0.05 is close to that employed in the the- 
oretical calculations used to generate Fig[4j the onset of 
the yield stress clipping effect, already manifest in Figj4j 
is not clear in the form of a(t) shown in Figj9j It would 
therefore seem likely that considerably smaller values of 
the Peclet number are required to observe this effect in 
our two-dimensional simulations. Nevertheless, the gen- 
eral form of the nonlinear stress is very similar, on a qual- 
itative level, to that predicted by the schematic model in 
Figj4] Both simulation and theory exhibit a flattened 
and asymmetric peak which is skewed to the left. 

In order to analyze more closely the stress signal we 
show in Figs lOj 11 and [12] the fundamental coefficients, 



higher harmonic intensities and phase shifts, respectively. 
The dependence of G' n and G'^ on 70 is strongly reminis- 
cent of that predicted by the schematic model (FigM. 
Within the range 0.01 < 70 < 0.06 the system begins to 
deform plastically leading to a reduction in G[ and an 
increase in G'{ (and cross at 70 = 0.05), reflecting the in- 
creasing importance of dissipative processes. The height 
of the peak in G'{ is rather less pronounced than that 
predicted by the schematic model. Beyond 70 = 0.06 
both Gi and G'{ decrease monotonically, exhibiting an 
asymptotic power law decay which is well described by 
the power laws G" ~ 7,7 - 68 and G[ ~ % 1A5 . Although 
these values do not satisfy perfectly the empirical rela- 
tion 1/ — 2v, the deviation of the simulation exponent 
ratio v' jv = 2.13 may well be attributable to numerical 
error. Despite this discrepancy, the absolute value of the 
exponents v = 0.68 and v' = 1.45 compare well with 
those emerging from the schematic model (y = 0.65 and 
v' = 1.3). The schematic model considered in the present 
work would thus seem to be more realistic than either a 
simple Maxwell model [v — 1 and v' = 2) or the micro- 
scopic MCT approach of Miyazaki et al. {v — 0.9 and 
v' = 1.8), at least on the basis of our simulation results. 

In Fig{TT] we show the intensities of the third and fifth 
harmonic as a function of 70- Higher order terms were 
found to be highly susceptable to the effects of statis- 
tical noise in the simulation data and have thus been 
omitted. Upon increasing the strain amplitude beyond 
70 = 0.03 the system leaves the linear response regime 
and the contribution of the third harmonic grows. For 
strains exceeding around 0.1 the fifth harmonic also be- 
gins to play a significant role in determining the stress 
response. In keeping with the schematic models predic- 
tions, both I3/I1 and I5/I1 exhibit a maximum, albeit 
more sharply peaked and shifted to slightly larger strain 
values approaching unity. The corresponding phase shifts 
also share the general features of the schematic model 
predictions, in particular, for large values of 70 we find 
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FIG. 13: (Experiment) Symbols: The experimentally mea- 
sured flowcurves for three different temperatures, T = 20 
{4> c « = 0.57), T = 18.4 {<j) cS = 0.60) and T = 15.1 
(0cff = 0.65). Lines: Theoretical fits to the data using the 
parameters in table [I] 
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FIG. 14: (Experiment) Symbols: The experimentally mea- 
sured linear response moduli for three different temperatures, 
T — 20 (cj> cH = 0.57), T = 18.0 {<f> eS = 0.60) and T = 15.0 
{4>cs = 0.65). Lines: Theoretical fits to the data using the 
parameters in table [I] 



that 61,63 and £5 saturate to ir/2, 3ir/2 and 5ir/2, respec- 
tively. These results suggest that the series ( 13 1 reduces 
to 

OO 

=7o J] WiH(-l) n cos((2n+l) W t), (27) 

n=l 

for large values of 70 . 



C. Experimental results 

1. Flow curves and linear response moduli 
In Fig. [13] we show the experimentally measured flow 



curves and in Fig. 14 the linear response moduli for three 
different temperatures (corresponding to three different 
volume fractions). Reduced units are employed for both 
the control parameters 



Pen 



D ' 



Pe = ^ 
" Do ' 



(where D is obtained from the solvent viscosity using 
Stokes' law) as well as for the shear stress and the moduli 
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For small Peclet numbers the flow curve measured at 
20° C (corresponding to cf> c g = 0.57) shows a first Newto- 
nian plateau. As Peo is increased we observe a decrease 
in viscosity, followed by a second Newtonian plateau, 
both of which are typical for a shear thinning fluid. The 
viscoelastic character of the sample at 20° C is demon- 
strated by the linear response moduli shown in Fig 14 



For intermediate frequencies G' and G" cross, indicating 
the presence of a structural a-relaxation process. 

The flow curve measured at the lower temperature 
18° C (4> c ff — 0.60) displays a more pronounced plateau 
region. However, for the lowest frequencies investigated 
a slight decrease from the plateau is evident, suggesting 
the existence of an a relaxation time which has shifted 
out of the experimental frequency window. The corre- 



sponding linear moduli shown in Fig 14 show a distinctive 
plateau region of G' at intermediate values of Pe u , fol- 
lowed by a decrease for small Pe u values, consistent with 
the existence of a crossover point and, therefore, a fluid 
relaxation. This is expecially apparent in G"(u>) which 
continues to rise as the frequency is decreased. Extrapo- 
lation of the measured data to lower frequencies suggests 



Pe" 



10" 6 - 10~ 7 , where G' and G" cross. For the 



lowest temperature investigated, 15° C (0 e ff = 0.65), the 
flow curve exhibits a constant plateau down to the lowest 
values of Pe and the storage moduli remains constant at 
low Pe^. The sample at 15° C may thus be considered 
as a glass for which additional 'hopping' processes lead 
to an increase in G" at low frequencies. 

The procedure by which experimental data may be 
fit using the schematic Fj 2 model of [57] is already well 
documented [32]. Fitting the experimental data for the 
flow curve and linear moduli using the present schematic 
model proceeds analogously. For a given volume fraction, 
a fixed set of model parameters may be found which fit 
both the flow curve and the linear moduli. It is thus 
possible to determine the separation parameter e, the 
vertex v a and the decay rate r. The parameter j c is ob- 
tained as an additional parameter for the description of 
the flow curve. The high frequency viscosity rf^ is only 
important for the frequency spectrum (and is connected 
with the high shear viscosity vj^ via 77^ = 77^ + v a j (2.T) 
[31)). Fixing the parameters by these two experiments 
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FIG. 15: (Experiment) Strain sweeps at three different fre- 
quencies for the temperature 15.0°C. Circles: Experimen- 
tal data. Lines: Theoretical fits using the parameter set 
(e, Vg, r, 7e, x) obtained by simultaneously fitting the flow 
curves (Fig |13[ ) and linear response moduli (Fig |14[ ). 



in the linear viscoelastic and the steady state determines 
all information needed to calculate the nonlinear oscil- 
latory behaviour (parameters are summarized in table 
|l|). Therefore, the experimental data sets of the defor- 



mation sweeps (Figs 
tests (Figs 



18 



15 - 17 1 and the oscillatory time 



and 19) are solely described by the the- 



II and the fixed parameter sets obtained 



14 and 13 without further modification. 



ory of Section 
by fitting Figs 
In this sense, the theoretical results to be presented for 
large amplitude oscillatory shear are predictions, as no 
further fitting is required. 



FIG. 16: (Experiment) Strain sweeps at three different fre- 
quencies for the temperature 18.0° C. Circles: Experimental 
data. Lines: Theoretical predictions. 



employing the parameter sets determined by the fitting 
procedure described above. 

In Fig |15| we show the results of strain sweep experi- 
ments at 15° C (4> c ff — 0.65) for three different values of 
Pe^. For strains up to around 1% the system shows lin- 
ear response behaviour, beyond which dissipation starts 
to increase, leading to a growth of G" up to a maximum 
in the range 10-20% strain. The growing dissipation and 
eventual crossing of G' and G" as a function of 70 indi- 
cate the breaking of microscopic particle cages. 

For higher deformations G' and G" display a power 
law decay. The exponents v and v' obtained at different 
temperatures and frequencies are given in table In] 



2. Nonlinear response 

The nonlinear regime has been tested using two dif- 
ferent experiments: deformation sweeps at 1, 0.1 and 
0.01 Hz (see Figs. 15 - 17 ) and oscillatory time series mea- 
sured at 1, 0.1 and 0.01 Hz for various deformation am- 
plitudes ranging from the linear to the non-linear regime 
(see Figs. 18 and 19). The complete set of nonlinear os- 



T [°C] 


0eff 


e 
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7c 




20.0 


0.57 


-2.45 x 10 -3 


59 


100 


0.18 


42 


18.0 


0.60 


-2.2 x 10~ 4 


85 


100 


0.19 


36 


15.0 


0.65 


5 x 10" 5 


115 


105 


0.28 


24 



cillatory data is solely described by the schematic MCT 



TABLE I: Schematic model parameters obtained by fitting 
experimental data for the flow curves (Fig 13 I and linear re- 
sponse moduli (Fig |14[ ). These parameters are then used to 
make theoretical prediction for large strain amplitude oscilla- 
tory experiments. 
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FIG. 17: (Experiment) Strain sweeps at a fixed frequency 
Pe„ = 0.01917 for the temperature 20°C. Circles: Experi- 
mental data. Lines: Theoretical predictions. 



The theoretical predictions shown in Fig jl5| are in ex- 
cellent agreement with the experimental data and cap- 
ture both the height and location of the maximum in G' 
rather well, although the departure from linear response 
appears to be less abrupt than in experiment, indicating 
a more gradual breaking of cages with increasing 70 . We 
note that the discrepancy between theory and experiment 
in the linear response regime for the lowest frequency 
considered has its origins in the linear moduli fits pre- 



sented in Fig 14 For glassy states there occur additional 
physical relaxation processes at low frequency which are 
manifest in an upturn of the linear response G" {uj) at low 
frequencies and which are not captured by mode-coupling 
based theoretical approaches. Particularly significant is 
the excellent agreement between theory and experiment 
in the large strain regime for which the cage structure 
has been broken up by the applied flow. The power law 
decay of the experimental data is well described by the 
theoretical exponents v = 0.65 and v' = 1.3. 

Both the numerically obtained theoretical results and 
experimental data shown here demonstrate the exponent 
relation v' « 2u (as in the theory of [40]). It is interesting 
to note that although this relation does not appear to be 
truly universal, broadly similar behaviour has been found 
for a variety of different materials, all of which show the 
deformation behaviour classified by Hyun et al. [59j as 
type III (weak strain overshoot). A few examples are 
e.g. a Xanthan gum solution [7] with v' = —1.53, v = 0.64 
and a ratio v' jv — 2.4, anchor spreadable butter and 
promise spread which yield exponents —2.1 < v' < —2.0 
with v = —0.9 and a ratio of v' jv = 2.3, or the hard- 
sphere solution of Miyazaki et al. 40 (PMMA spheres of 
197nm in a mixture of decaline and cyclohcptylcbromide) 
which show v = 0.7 and v' — 1.4. 

In Figs |16| and [17] we show further strain sweep mea- 
surements for temperatures 18° C (</> e ff = 0.60) and 20° C 
(0cff = 0.57). The overall level of agreement between 
theory and experiment is even better than that found 
for T=15° C. In particular, the strain sweep measured 
at 20° C is very well described by the theory and, taken 
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FIG. 18: (Experiment) The stress response measured in 
LAOS experiments for strain amplitudes from 70 = 0.03 to 
70 = 5 (black line) and the associated Lissajous figures illus- 
trating the nonlinear character of the response. The exper- 
iments are performed at T = 15.1°C a glassy state and at 
a frequency of 1Hz (corresponding to Pe^ = 0.02533). At 
70 = 0.03 the response is almost entirely elastic, emphasising 
the proximity of the quiescent state to the glass transition. 
At 70 = 5 the system is almost purely viscous. The increase 
in dissipation with increasing 70 is reflected in the increasing 
area enclosed by the closed Lissajous curves. The yield stress 
is indicated by the broken red lines. Theoretical results are 
given by broken blue lines. 



together with the results shown in Figs 13 and 14 demon- 
strate the accuracy with which the present schematic 
model may be used to describe the flow curve, linear 
moduli and strain sweep data of a dense colloidal fluid 
using a single fixed set of model parameters. 

For the oscillatory time series the schematic model cal- 
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FIG. 19: (Experiment) As in Fig 18 but at a frequency of 
0.1Hz (corresponding to Pe^ = 0.002533). Full black lines: 
Experimental data. Broken blue lines: Theoretical results. 
The yield stress is indicated by the broken red lines. 



culations and the simulation are in good agreement with 
the experimental data. The direct comparison of theory 
and experiment for two frequencies (1 Hz and 0.01 Hz) re- 
spectively the Peclet numbers (Pe w = 0.025 and 0.0025) 
is given in Figs. [18] and [T9] For small deformations, a 
linear viscoelastic behaviour is indicated by the nearly 
perfect sinusoidal output signal, but becomes distorted 
as the strain amplitude is increased. For Pe u = 0.025 
the stress signal displays flattened asymmetric peaks at 
intermediate values of 70, consistent with a regime of cage 
breaking around 70 = 7c- In contrast to the MCT predic- 
tions, the data show a pronounced dip at the top of the 
asymmetric peak. For high deformations, the peak shape 
approaches a semi-spherical shape with a vanishing but 
still visible dip or overshoot at the beginning edge of the 



peak. For the smaller frequency at Pe^ — 0.0025 indi- 
cations of the effect of the merging a max and the yield 
stress are found. The peaks show in the intermediate and 
high deformation range a more cut off shape, although 
the dip still remains, in contradiction to the MCT. The 
peak shape for middle 70 drops faster in the experiment 
as for the more box-like shape of the MCT time signals. 
However this shape is found in the experiment for high 
70 as well. Alltogether the experiments and the theory 
fit well for all compared Peclet numbers and amplitudes, 
although some small deviations of the shape exist. The 
predictive character of the schematic MCT model for the 
oscillatory time series is remarkable, due to the fact that 
the shapes are not fitted but calculated with the param- 
eter set defined from the flow curves and the frequency 
test in the linear viscoelastic regime. 



3. Fourier 



In order to provide a more detailed analysis of the time 



series shown in Figs 18 and 19 we have calculated the in- 
tensities and phase shifts (see Eq.(13l). The phase shifts 
are calculated from the real and imaginary parts of the 
Fourier transformed signal. Representative FFT-spectra 



are shown in Fig. 20 for the measurements made at at 
15.1° C, 1Hz and for amplitudes 70 = 0.03 and 1. The 
spectrum demonstrates that at this low strain amplitude 
the higher harmonics do not contribute significantly and 
have an intensity only around 0.1% of the basic harmonic 
I\. Although peaks for even and odd harmonics are 
visible, their low intensities are not far away from the 
noise level. We therefore consider this measurement to 
be within the linear viscoelastic range. 

At the same conditions the MCT calculations show a 
rather different picture (see Fig. 20 d). Although the time 



signals for this 70 do not deviate strongly from the exper- 
imental equivalent (see Fig. 20 1 , the effect in the Fourier 
spectrum is noticable. Up to the 11th harmonic all odd 
harmonics can be clearly separated from the base line, 
with /3//1 around 1%. Furthermore, in the theoretical 
FFT-spcctrum an exponential decay of the intensities is 
apparent, a feature which first becomes evident in the 
exterimental data for 70 > 0.158. At the larger strain 
amplitude, 70 = 1, the MCT calculations and the ex- 
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0.65 
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-0.633 


-1.262 


1.99 


15.0 


0.65 
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-0.742 


-1.404 
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15.0 


0.65 


0.01 


-0.788 


-1.555 
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18.0 


0.60 
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-0.631 


-1.297 


2.06 


18.0 


0.60 


0.1 


-0.712 


-1.424 


2.00 


18.0 


0.60 


0.01 


-0.824 


-1.630 


1.98 


20.0 


0.57 


1 


-0.630 


-1.250 


1.98 



TABLE II: Experimentally measured values of the exponents 
v and v 1 dictating the decay of G' and G" as a function of 70 
for large values of 70. 
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FIG. 20: (Experiment and Theory) FFT-spectra of the oscil- 
latory time test signals. The experiment was performed at 
15.1° and 1Hz. a) (experiment) and b) (theory) are the FFT 
spectra at 70 = 0.03 and c) (experiment) and d) (theory) at 
an amplitude of 70 = 1. 



perimental FFT-spectra are very similar (see Fig. 20: 
and d). The present theory thus provides a qualita- 
tive description of the higher harmonics which becomes 
semi-quantitative as the strain amplitude becomes sig- 
nificant (70 ~ 1). The theory provides a sensible int- 
perpolation between linear response and large amplitude 
regimes, both of which are captured accurately. 

Some insight into the origin of the discrepancies at 
small excitation amplitude may be obtained by consid- 



ering the theoretical predictions for the buildup of stress 
following the onset of steady shear flow (see also [35]). 
Within the elastic regime (i.e. well below the 'stress 
overshoot' identifying the yield strain) the system should 
display Hookian behaviour with a clearly defined elastic 
constant. However, schematic model calculations of the 
stress for this protocol (either using the present model or 
the F^ 2 model of [37]) exhibit devaiations from Hookian 
behaviour at small strains. This feature of the schematic 
models is related to the slowness of the j3 decay of the 
transient density correlator onto the plateau. We can 
thus speculate that the discrepancy between the Fourier 
spectra of theory and experiment at low values of 70 
may be due to the excessively slow decay of the corre- 
lator to its plateau value, which is an inherent feature 
of any model based on the original schematic F 12 model 



In Fig. [21ft the amplitude dependence of the exper- 
imentally measured third harmonic is shown for three 
different frequencies in the glass (15.1° C). It can be seen 
that as the frequency is reduced, the onset of the non- 
linear regime, indicated by increase of the normalized 
third harmonic intensity, moves to lower values of 70. In 
Fig. |21b we show the same quantity at a fixed frequency 
of 1Hz for three different temperatures. Surprisingly no 
strong influence of a change in the volume fraction is 
found, apart from a small deviation in the onset of the 
non-linear regime in the glassy state, which is shifted to 
higher deformations. The only significant deviation for 
the volume fractions considered occurs at intermediate 
deformations, as the starting and end values coincide. 

The phase shifts are given for the measurement at 1 Hz 
(Pe^ = 0.025) and 15° C in Fig. [2lk In the case of the 



fundamental Si the phase shift is found as expected to 
start at 0° for small amplitudes and to end at 90° for 
high deformations, which corresponds to a cosine. The 
curve progression of the experimental data and the MCT 
calculation fits perfectly for Si. The theory deviates for 
the phase shifts of the higher harmonics due to the exces- 
sive contribution of non-linearities at lower 79. However, 
the limiting values for high deformations are found to 
coincide (f for S%, 4^ for S3 and ^ for <$ 5 ). 

We have also included the data from our two- 
dimensional simulations into Figs. [21] : and d. In order 
to obtain a reasonable comparison we found it necces- 
sary to empirically multiply the strain employed in the 
simulations be a factor of three, 70 = 3 • 7q 4 "\ That 
such an empirical rescaling is necessary is not surprising, 
given that only qualitative comparison is to be expected 
when comparing three dimensional experimental results 
with those of two dimensional simulations. It is therefore 
gratifying that the phases of the rescaled ^ zm are found 
to describe the experimental data very well, not only for 
£1 but also for S3 and S5. 

The experimentally measured I3/I1 at 15.1° C and 
1 Hz is found to indicate the onset of the transition from 
linear to non-linear regime at around 70 = 0.04 — 0.05. 
This can be seen by comparing with the strain sweep 
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FIG. 21: (Experiment and Theory) a) I 3 /h at 15.1° C for 
different frequencies (red circle: 1Hz, black square: 0.1 Hz 
and blue triangles up: 0.01 Hz), b) Temperature dependence 
ofh/h at 1Hz: 15.1° C (blue squares), 18.4° C (red circles) 
and 20.9° C (green triangles down), c) Experimental phase 
shifts (symbols) for 15.1°C at 1Hz. The simulation data is 
plotted versus 70 = 3~(o m as thin lines. The MCT calcu- 
lations is given in thick lines. The phase shifts are Si (red 
triangles up and dotted lines), 63 (blue circles and dashed 
lines) and S$ (black squares and solid lines), d) Normalized 
intensities at 15.1° C and 1Hz; symbols mark the experimen- 
tal data, thick lines the FFT of the MCT time signals and 
thin lines the shifted simulation results: I2/I1 (green trian- 
gles down and solid line), I3/I1 (red triangles up and dotted 
lines), I5/I1 (blue circles and dashed lines) and I7/I1 (black 
squares and solid line). 



data. It shows the maximum of G" at 70 ~ 0.16, which 
is connected with the breaking of the cages. This value 
could be correlated with the raising of above the 

'noise' level of 0.1 • The description of I3/I1 is pos- 

sible with the formula obtained by Wilhelm (13/1(70) = 
A [l - (1 + (-B7o) c ) -1 ]) 0. For all experimental data 
the exponent c or the slope of the increase is found to be 
in the range of 2.2 — 2.7. The rescaled simulation results 
for the harmonics describe the experiment rather closely, 
whereas the MCT calcuations are found to show a dif- 
ferent behaviour. All odd harmonics start at low 70 at 
higher values, whereas for high deformations the inten- 
sities coincide with the experimental results. Thus, also 
the slope of I3/I1 in this diagram changed to c ps 1.7 
instead of the experimental c = 2.7. These discrepancies 
between experiment and theory may well be attributable 
to the slow decay of the correlator to the plateau, as dis- 
cussed above. 

Whereas the mode-coupling calcuations show a very 
early transition from the linear to the nonlinear regime 
at 7 w 0.01 (see Figs. 15 to 17 1, for reasons discussed 



above, the experiments for the frequencies and tempera- 
tures investigated show a deviation from linear behaviour 
at 7 ps 0.04 — 0.05. For the rescaled simulation results 
of Fig. [lOj the deviation starts at 7 ps 0.03. The on- 
set of the yielding process is correlated with the onset 
of higher harmonics, starting with the third harmonic 
(see Fig. 21 i). Although the difference of the onset and 



intensity of the higher harmonics between theory and ex- 
periment is significant, the time signals are only slightly 
influenced (Fig. 18 and 19). The FFT is a very sensitive 
method to analyse the time signals, so very small devia- 
tions in the time signal can cause a remarkable difference 
in the spectra. Increasing the strain amplitude results in 
more asymmetric time signals with maxima and minima 
shifted to the left. 

It was found in experimental, theoretical and simu- 
lation strain sweeps that the maximum of G'[ lies very 
close to the crossing point of G[(^o) and G'{(-fo). For 
the glassy sample {(f) = 0.65) this point was located at 
7 ss 0.15. Beyond this value, the fifth harmonic be- 
gins to increase, which can be seen in Fig. [21p l for the 
experimental and rescaled simulation data. The theo- 
retical time signals do not show such a sharp transition 
from the linear to the nonlinear regime. Furthermore, 
the simulation data starts at higher intensity ratios as 
the experiment; an effect even more pronounced for the 
theoretical data. The schematic model thus predicts a 
more gradual transition between solid and fluid. The 
onset of the fifth harmonic was also found for microme- 
ter sized particles [HI] to be correlated with the crossing 
point of G[ and G'{. The view that the yield strain is 
indicated by the onset of the fifth harmonic is supported 
by the results of Petekidis [B2]. For larger deformations, 
the time signals show a strain softening behaviour, which 
is typical for shear thinning fluids. The phase shifts of 
experiment, simulation and theory approach the limiting 
values of n -tt/2, where n is the index of the harmonic. In 
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addition, the harmonics for theory and experiment show 
the same limiting values lat large 70 for each harmonic 
(21% for I z /h, 10% for I 5 /h, 6% for hjh and 4% for 
/9/ii). Furthermore the theoretical harmonics exhibit a 
slight decrease at the highest calculated strains, which is 
not observed in the experimental data, perhaps due the 
choice of a too small measurment range. In the simu- 
lation the decrease of the intensities is much more pro- 
nounced for high deformations, as the intensities do not 
show a plateau but only a maximum (at 18% for J3//1 
and at 8% for I 5 /h). 



VII. CONCLUSION AND OUTLOOK 

We have used a combination of theory, experiment and 
simulation to investigate the nonlinear stress response of 
dense colloidal dispersions under large amplitude oscil- 
latory shear flow. The theory employed is a recent ex- 
tension [23] of the well studied F7 2 model [27] ■ A key 
physical mechanism captured by the theory is the yield- 
ing of local particle cages and the subsequent onset of 
nonlinearity. In contrast to a related approach presented 
by Miyazaki et al. [35] HO]; the theory employed here 
makes predictions for the higher harmonic contributions 
to the stress signal and thus enables the extent of the 
linear response regime to be addressed. In order to make 
contact with experiment, the theory requires only the 
simultaneously determined parameters from the steady 
state flow curve and the frequency dependence of the 
linear moduli to make predictions for the nonlinear vis- 
coelastic response. We have thus compared the theory 
with rheological experiments performed on concentrated 
suspensions of thermo-sensitive core-shell particles. 

The first nonlinear experiment considered was a strain 
sweep, performed at three different volume fractions and 
frequencies. Here we found very good agreement between 
experiment and theory, although there are small devia- 
tions concerning the onset of the nonlinear viscoelastic 
regime. The values of G\ and G'{ in the linear regime, the 
crossing point of G[ and G'{, maximum of G'{ and asymp- 
totic large strain behaviour of G[ and G'{ are all well 
described by the theory. Brownian simulation of a two- 
dimensional system of hard discs enable only a qualita- 
tive comparison, but show behaviour broadly consistent 
with our experimental data, indicating that the yielding 
process is not strongly dependent upon either the mate- 
rial details or dimensionality of the sample. Moreover, 
the simulations do not contain hydrodynamic interac- 
tions, suggesting that these are not of great importance 
for yielding. The ratio of the slopes of G' x and G'[ for 
large 70 in the shear molten state yields the exponent ra- 
tio v' jv = 2.0±0.1 in theory simulation and experiment. 
The fact that various other materials (of type III in [59] ) 
display similar asymptotic behaviour may suggest a uni- 
versal mechanism underlying the oscillatory response of 
shear molten viscoelastic fluids. 

The second type of experiment considered were os- 



cillatory time sweeps. In this case the cage yielding is 
expressed by the deviation of the signal from a sinu- 
soidal form, showing a characteristic asymmetric peak 
in the yielding regime, which is followed in the shear 
molten state by a typical strain softening semi-circular 
peak shape. The agreement of the time signals obtained 
from theory, simulation and experiments is good. How- 
ever, small deviations in the time series obtained using 
the three methods lead to stronger deviations in the pa- 
rameters of the Fourier transformed time signals. This 
serves to emphasize the fact that FT-rheology is a very 
sensitive method, sensitive on the logarithmic scale, ca- 
pable of detecting e.g. the fifth harmonic, to an accuracy 
of less than 1 promille. This analysis has shown that the 
onset of the third harmonic heralds the start of the non- 
linear regime and that the maximum in G", which here 
approximately coincides with the both crossing point of 
G' x and G'{ and the yield strain are correlated with the 
onset of the fifth harmonic. Moreover, the plateau val- 
ues of the phase shifts at high deformations are found to 
follow n ■ 7r/2 with n beeing the index of each harmonic. 

Despite the good overall level of agreement between 
theory and experiment there remain aspects which could 
be improved. Firstly, the extent (in 70) of the linear re- 
sponse regime is apparently too small within the present 
theory, with consequences for the variation of the higher 
harmonic intensities with amplitude. As noted, this fail- 
ing has its origins in the quiescent transient density corre- 
lators predicted by the F12 model, upon which our more 
recent schematic model is based. It would thus be de- 
sireable to improve this fundamental aspect of the the- 
ory. Secondly, the stress response measured in experi- 
ment displays a more asymmetric waveform than that 
predicted by the schematic model. This aspect can po- 
tentially be connected to the fact that our model, when 
applied to calculate the buildup of stress upon the onset 
of shear, does not generate a stress overshoot. While a 
stress overshoot does occur in approximate solutions of 
the fully microscopic mode-coupling expressions [35] , this 
aspect is lost in making the ansatz ( 16 1 for the shear mod- 



ulus, as the strain dependent vertex functions appear- 
ing in the microscopic theory are replaced by a constant 
V a . Empirically reincorporating a strain dependence into 
the theory via the replacement v a — > u CT (7o)) such that 
the overshoot is recovered, may also lead to an increased 
waveform asymmetry in the stress response and, thus, 
better agreement with experiment. However, when mod- 
ifying the schematic model in this fashion, care must be 
taken not to destroy its existing positive features. 

To summarize, it can be concluded that the 
schmematic mode-coupling model 24J can make accurate 
predictions in the nonlinear viscoelastic regime, based 
purely on parameters fixed by the steady state stress 
and linear viscoelastic behaviour. This constitutes the 
first truly time-dependent test (other than step strain) 
of the schematic model proposed in [24]. Fourier trans- 
form analysis of time series obtained for various strain 
amplitudes and frequencies provides a wealth of experi- 
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mental information regarding the mechanical response of 
a material. As our present experimental setup enables 
the investigation of transient flows, we anticipate that 
the study of such flow protocols, in combination with 
the present results, should enable a complete rheological 
characterization of our colloidal system. 
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